/**
 * 
 */
package br.unesp.ibb.lbbc.fba;

import java.io.IOException;
import java.util.ArrayList;
import java.util.Collection;

import javax.xml.stream.XMLStreamException;

import org.apache.commons.math3.optim.MaxIter;
import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.linear.LinearConstraint;
import org.apache.commons.math3.optim.linear.LinearConstraintSet;
import org.apache.commons.math3.optim.linear.LinearObjectiveFunction;
import org.apache.commons.math3.optim.linear.NonNegativeConstraint;
import org.apache.commons.math3.optim.linear.Relationship;
import org.apache.commons.math3.optim.linear.SimplexSolver;
import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;


import org.sbml.jsbml.ListOf;
import org.sbml.jsbml.LocalParameter;
import org.sbml.jsbml.Model;
import org.sbml.jsbml.Reaction;
import org.sbml.jsbml.SBMLDocument;
import org.sbml.jsbml.SBMLReader;
import org.sbml.jsbml.Species;
import org.sbml.jsbml.SpeciesReference;

import br.unesp.ibb.lbbc.util.FileEC;
import br.unesp.ibb.lbbc.util.MathTools;

/**
 * @author Esther
 * 
 */
public class Stoichometry {

	/**
	 * Get stoichometric matrix from SBML file
	 * 
	 * */
	public double[][] calculaStoichometry(String pathName) {

		Collection constraints = new ArrayList();

		java.io.File sbmlFile = new java.io.File(pathName);
		SBMLReader reader = new SBMLReader();
		SBMLDocument doc = null;
		try {
			doc = reader.readSBML(sbmlFile);
		} catch (XMLStreamException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}

		Model model = doc.getModel();
		ListOf<Reaction> listReactions = model.getListOfReactions();
		ListOf<Species> listMetabolites = model.getListOfSpecies();

			
		int reactionSize = listReactions.size();
		int metabSize = listMetabolites.size();

		double[][] matrix = new double[metabSize][reactionSize];

		for (Reaction reaction : listReactions) {

			double[] upperBound = new double[reactionSize];
			double[] lowerBound = new double[reactionSize];

			ListOf<SpeciesReference> reactants = reaction.getListOfReactants();
			for (SpeciesReference ref : reactants) {
				int a = listMetabolites.indexOf(ref.getSpeciesInstance());
				int b = listReactions.indexOf(reaction);
				matrix[a][b] = -1 * ref.getStoichiometry();
			}

			ListOf<SpeciesReference> products = reaction.getListOfProducts();
			for (SpeciesReference ref : products) {
				int a = listMetabolites.indexOf(ref.getSpeciesInstance());
				int b = listReactions.indexOf(reaction);
				matrix[a][b] = ref.getStoichiometry();
			}
			
			
			LocalParameter lowParam = null;
			
			try {
				lowParam = reaction.getKineticLaw().getLocalParameter("LOWER_BOUND");
			} catch (Exception e) {
				// TODO Auto-generated catch block
				//e.printStackTrace();
			}
						
			if (lowParam != null) {
				lowerBound[listReactions.indexOf(reaction)] = 1;
				constraints.add(new LinearConstraint(lowerBound,Relationship.GEQ, lowParam.getValue()));
			}
			//else {constraints.add(new LinearConstraint(lowerBound,Relationship.GEQ, -10000.0));}
			
			
			LocalParameter upperParam = null;
			
			try {
				upperParam = reaction.getKineticLaw().getLocalParameter("UPPER_BOUND");
			} catch (Exception e) {
				// TODO Auto-generated catch block
				//e.printStackTrace();
			}
			
			if (upperParam != null) {
				upperBound[listReactions.indexOf(reaction)] = 1;
				constraints.add(new LinearConstraint(upperBound,Relationship.LEQ, upperParam.getValue()));
			}
		//	else {constraints.add(new LinearConstraint(upperBound,Relationship.LEQ, 10000.0));}
		}

		
		//MATRIZ ESTEQUIOMETRICA
		for (int i = 0; i < matrix.length; i++) {
			constraints.add(new LinearConstraint(matrix[i], Relationship.EQ, 0));
		}
		
		//FUNCAO DE BIOMASSA 
		Reaction biomass = model.getReaction("R_BiomassEcoli");
		//Reaction biomass = model.getReaction("R_Ec_biomass_iJO1366_WT_53p95M");
		
	//	FileEC.saveFile(matrix, "matrixJeni");
		
		double[] max = new double[reactionSize];
		
		max[listReactions.indexOf(biomass)]=1;
	
		
		
		//double [] max = {0,0,0,1,1,1};
		
		LinearObjectiveFunction f = new LinearObjectiveFunction(max, 0);
		SimplexSolver solver = new SimplexSolver();
		
	    PointValuePair solution = solver.optimize(new MaxIter(1000),f, new LinearConstraintSet(constraints),GoalType.MAXIMIZE, new NonNegativeConstraint(false));		

		// get the solution

		for (int i = 0; i < solution.getPoint().length; i++) {
			System.out.println(solution.getPoint()[i]);
		}

		double min = solution.getValue();
		System.out.println("\n" + min + "\n");

		return matrix;
	}

	public static void main(String[] args) {
		Stoichometry s = new Stoichometry();
		//String pathName = System.getProperty("user.dir") + "\\ecoli.xml";
		String pathName = System.getProperty("user.dir") + "\\SBML_Jenniffer904.xml";
	//	String pathName = System.getProperty("user.dir") + "\\small.xml";
		double[][] matrix = s.calculaStoichometry(pathName);	
	}

}
